Tarea 3
Tarea 3: Frequentist Inference II
Integrantes :
- Nicolás García Ríos
Cuerpo Docente:
- Profesor: Felipe Bravo M.
- Auxiliar: Sebastian Bustos e Ignacio Meza D.
Fecha límite de entrega:
Índice:
Objetivo
Bienvenid@s a la tercera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en el diseño de experimentos, test de hipótesis y regresión lineal. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Instrucciones:
- La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacer de forma individual.
- La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
- El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
- Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
- No serán revisadas tareas desarrolladas en Python.
- Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
- Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.
Primera Parte: Preguntas Teóricas
A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas de forma breve.
Pregunta 1:
Determine si las siguientes regresiones son lineales para los parámetros \(\beta_{i}\)
\(y(x) = \beta_{0} + \beta_{1} x\)
\(y(x) = \beta_{0} + \beta_{1} x + (\beta_{3}+x)^{3}\)
\(y(x) = \ln(x^{2}+\beta_{0})\)
a): Tenemos un \(\beta_{0}\) sin ninguna variable independiente y \(\beta_{1}\) acompañado de \(x\). Luego como el grado de la ecuación es 1, esta regresión es lineal.
b): Expandiendo la última expresión, tenemos que \[y(x) = \beta_{0} + \beta_{1} x + (\beta_{3}+x)^{3}\] \[y(x) = \beta_{0} + \beta_{1} x + \beta_{3}^{3}+3\beta_{3}^{2}x + 3\beta_{3}x^{2} + x^{3}\] \[y(x) = (\beta_{0} + \beta_{3}^{3}) + (\beta_{1} + 3\beta_{3}^{2}) x + 3\beta_{3}x^{2} + x^{3}\]
Pero como la ecuación es de grado cúbico, tenemos que esta es una regresión polinomial de grado 3, no lineal.
c): Esta no es lineal, ya que sus parámetros están dentro de una función logarítmica y por ende no lineal.
Pregunta 2:
Una universidad esta interesada en saber cuantos alumnos gustan del anime, para esto realiza una recopilación de datos y llega a la construcción del siguiente modelo lineal simple:
\[\hat{N°\_de\_fanaticos\_del\_anime}=300+90*numero\_de\_semestres\]
¿Como podemos interpretar el intercepto y pendiente del modelo?, ¿El intercepto y la pendiente tienen una interpretación coherente para cualquier modelo lineal?, si no es así, de un ejemplo.
Podemos pensar, según este modelo lineal, que en 0 semestres tenemos 300 fanáticos de anime, y por cada semestre llegan 90 nuevos fanáticos. Esto por su puesto no tiene ninguna interpretación coherente, ya que es esperable que en el semestre 0 no exista ningún alumno matriculado en la universidad.
Con esto tenemos un ejemplo de que el itercepto y la pendiente de un modelo lineal no siempre tienen una interpretación u implicancia lógica. Otro ejemplo de esta discordancia es con el dataset visto en clases (datos antropológicos de una tribu), en ese caso vimos un modelo de regresión lineal simple con la altura dependiente de la edad, donde se obtuvo que un humano de 0 años mediría según el modelo unos 113 centímetros, y que cada año se estimaba un aumento de 0.9 centimetros en la altura. No es lógico tener que un recién nacido mida más de un metro, o que un anciano de 100 años mida más de dos metros.
Claramente estos tipos de modelos están pensados para que funcionen con datos promedios, o con datos generados explícitamente con funciones lineales. Aquí entramos nuevamente en la diferencia entre significancia estadística y real, dado que estadísticamente es relevante saber el intercepto y pendiente de un modelo lineal para poder ver el comportamiento de un grupo de datos, pero en la realidad estos datos no siempre son coherentes o importantes en la práctica.
Pregunta 3:
Considere un test con dos muestras no aparejado, explique porque se hace una corrección a los grados de libertad en el test de Welsh.
Se una corrección con el fin de incorporar al problema la diferencia entre los tamaños y las varianzas de las muestras, con el fin de penalizar estas diferencias. Así tenemos que con menor grados de libertad se tiene una mayor incertidumbre.
Pregunta 4:
Al realizar una regresión lineal simple con una variable categórica \((\beta_{0} + \beta_{1} \cdot \text{categorica})\) ¿que interpretación puede obtenerse del coeficiente que acompaña a la variable categórica?
Estamos trabajando cun una variable dummy, lo cual para variables categóricas tendremos que \(\beta_{0}\) es la media para el caso \(\text{categorica}=0\), lo que implica que \(\beta_{1}\) es la correción del modelo para \(\text{categorica}=1\), tomando como base a la media para \(\text{categorica}=0\). En otras palabras, la constante es el promedio para una de las variables categoricas, y la pendiente la diferencia entre las medias de estas.
Pregunta 5:
Discuta la siguiente frase:
" Hacer una regresión lineal mediante máxima verosimilitud requiere tener ciertas hipótesis probabilisticas de los datos, mientras que una regresión realizada mediante mínimos cuadrados no necesita tener ninguna hipótesis probabilista. "
Trabajar con un modelo basado en mínimos cuadrados solo necesita optimizar una función de predicción y reducir los errores, pero para trabajar con una regresión lineal mediante máxima verosimilitud es necesario tener una hipótesis o tener conocimiento respecto al comportamiento probabilístico de los datos. En otras palabras, utilizar máxima verosimilitud exige un conocimiento previo del comportamiento, mientras que mínimos cuadrados no.
Pregunta 6:
Explique porque el test de significancia sobre los parámetros de una regresión lineal se realiza bajo la hipótesis nula \(\beta_{H_{0}}=0\).
Esto es así debido a que este mide el grado de independencia de los parámetros de la regresión lineal (cuando \(\beta_{H_{0}}=0\), estos son independientes entre sí), así, lo que se busca es encontrar evidencia que permita descartar la hipótesis nula que asume independencia entre los parámetros y poder llegar a una relación lineal.
Pregunta 7:
¿Que nos dice la equivalencia de la máxima verosimilitud sobre los parámetros que componen una regresión lineal? ¿Qué nos permitirían calcular?
Nos dice que utilizar estimación de máxima verósimilitud sobre los parámetros de una regresión lineal es en realidad lo mismo que utilizar el modelo de mínimos cuadrádos, ya que parten de las mismas suposiciones iniciales (como comportamiento normal, u homocedasticidad), y llegan a la misma optimización. Esto nos permite calcular un modelo de regresión lineal ajustada en base a errores cuadráticos.
Pregunta 8
Consideremos una regresión lineal de una variable. En vez de mínimos cuadrados es posible minimizar la expresión
\[ \displaystyle{\sum_{i=1}^{n}}|y_{i}-\beta_{0}-\beta_{1}x_{i}| \] Explique en que se diferencia con mínimos cuadrados, de una posible ventaja y desventaja de este método (en comparación a mínimos cuadrados).
Esta se diferencia con SSE en el factor cuadrático del error, por lo que esta optimización es menos estricta con los errores a diferencia de los minimos cuadrádos que “penaliza” de mayor manera a los valores con mayor error, al aumentar cuadráticamente su valor.
Una ventaja de este modelo puede ser su computabilidad y robustez ante outliers. A mayor cantidad de datos el modelo cuadrático podría demorarse mucho más al tener que calcular el cuadrado de cada error en vez de trabajar con ellos directamente, y a demás, al no penalizar valores extremos evitamos que estos afecten gravemente al modelo.
Por otro lado, una desventaja podría ser su menor restricción. Utilizar SSE hace que el modelo tenga que minimizar la cantidad de “errores grandes”, por lo que es dificil tener grandes margenes de errores en un modelo de regresión. Esto por su puesto no ocurre con este nuevo modelo, ya que en vez de minimizar la cantidad de errores grandes, este minimizaría los errores sin importar que tan distintos sean del valor original, lo cuál afectaría al modelo si es que este no contiene outliers.
Pregunta 9:
Explique porque el coeficiente \(R^2\) tiende a crecer con el numero de variables.
Esto sucede porque a medida que aumenta la cantidad de variables en el conjunto de datos, el modelo reducirá naturalmente el error cuadrático con el fin de ajustar lo más posible el modelo con los datos, lo que se traduce en un aumento de la variable R2.
Pregunta 10
Un estudio de cáncer realizado por la institución X ha señalado que las personas que beben café poseen mayores probabilidades de padecer algún cáncer pulmonar. El estudio causo un gran revuelo en la población, por lo que una segunda institución ha decidido replicar el experimento, llegando a la conclusión que las personas que toman café tienden a fumar cigarrillos mientras beben esta bebida. Señale que tipo de variable serían los fumadores de cigarrillos en el estudio de cáncer pulmonar y explique cual es la característica de estas variables.
Los fumadores de cigarrillos son claramente una variable de confusión. Estas variables hacen que dos variables mutuamente relacionadas entre sí sean confundidas por causalidad en vez de correlación, y debido a que estas están muy ligadas entre ellas, se hace difícil observar la relación real.
Segunda Parte: Elaboración de Código
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:
# Manipulación de estructuras
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.1.0 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(dplyr)
library(tidyr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout# Manipulación de varios plots en una imagen.
library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.0.5##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(grid)Z-test
En clases se han visto diferentes tipos de test de hipótesis para demostrar una proposición sobre algún parámetro. Uno de los test vistos en clases es el Z-Test, el cual su distribución del test estadístico bajo la hipótesis nula se puede aproximar a una Gaussina. Para la aplicación de este test, resaltan los siguientes puntos:
- Cada uno de los puntos de la muestra deben ser independientes unos de otros.
- Al utilizar una distribución normal en la hipótesis nula, este test debería utilizarse cuando se tiene un número considerable de observaciones, ya que la sampling distribution de la media tiende a una gaussiana, de lo contrario se debería usar un T-test.
Para calcular la significancia estadística al igual que con otros métodos esta se debe calcular como:
- Menor/Cola-Izquierda (one-tailed): La Hipótesis Nula H0: \(\mu \geq \mu0\) vs Hipótesis Alternativa H1: \(\mu < \mu0\).
- Superior/Cola-Derecha (one-tailed): La Hipótesis Nula H0: \(\mu \leq \mu0\) vs Hipótesis Alternativa H1: \(\mu > \mu0\).
- Dos-Colas/Two-tailed: Hipótesis Nula H0: \(\mu = \mu0\) vs Hipótesis Alternativa H1: \(\mu \neq \mu0\).
Luego, dependiendo del objetivo del test tenemos las metodologías one-sample y two-sample. Utilizaremos One-Sample cuando nuestro objetivo es comparar la media de una muestra con la media de la población. El Z-score del One-Sample se define como:
\[Z-score_{One-Sample} = \dfrac{\bar x - \mu}{\dfrac{\sigma}{\sqrt n}}\] Donde \(\bar x\) es la media de la muestra, \(\mu\) es la media de la población, \(\sigma\) es la desviación estándar de la población y \(n\) es el tamaño de la muestra.
Por otro lado, se utiliza Two-Sample cuando queremos comparar la media de dos muestras. El Z-score de Two-Sample se define con la ecuación:
\[Z-score_{Two-Sample} = \dfrac{(\bar x_1 - \bar x_2) - (\mu_1 - \mu_2)}{\sqrt{\dfrac{\sigma_1^{2}}{n_1} + \dfrac{\sigma_2^{2}}{ n_2}} }\]
Donde \((\bar x_2 - \bar x_1)\) es la diferencia de las medias de la muestra, \((\mu_1 - \mu_2)\) la diferencia de las medias de la población, \(\sigma_{1,2}\) la desviación estándar de la población y \(n_{1,2}\) el tamaño de las muestras.
Multiples Test
En la práctica aparece la necesidad de testear múltiples hipótesis (por ejemplo en biología se pueden utilizar múltiples grupos de control o querer estudiar múltiples resultados de un mismo experimento), de esta forma la primera idea es testear individualmente cada una de las hipótesis, el problema de este enfoque es que la probabilidad de que se obtenga al menos un resultado significante crece rápidamente (con un nivel de significancia \(\alpha = 0.05\) y \(20\) test ya se alcanza una probabilidad de \(64\%\) de tener resultados significantes por azar).
Una forma de corregir los inconvenientes del método anterior es utilizar el método de Bonferroni correction quien propone cambiar \(\alpha\) por \(\alpha/m\) (donde \(m\) es la cantidad de test de hipotesis realizados), esto resulta que las probabilidades de rechazar por error se mantengan bajas. De esta forma los p-valores obtenidos en un test de hipótesis y al utilizar Bonferroni correction, quedan dados por el producto de un \(p-valor_{i}\) y la cantidad de test realizados: \(\text{p-valor}_{i}*m\).
Pregunta 1: “I´ve Got The Power!”
El objetivo de esta pregunta es programar la potencia de un test de hipótesis y observar como se comportan las la hipótesis nula v/s la alternativa para un Z-test. Con el desarrollo de este ejercicio, podrán visualizar las diferentes partes que conforman a un test de hipótesis, identificar que es el p-valor y evidenciar como varia la potencia de un test one-sample y two-sample al variar \(\alpha\).
Para recordar; sabemos que en estadística el concepto de potencia viene dado por:
\[Power = 1 - \beta\]
Donde \(\beta\) es la probabilidad de obtener un error de tipo II. Con esto, la potencia estadística viene a representar la probabilidad de rechazar la hipótesis nula cuando esta es falsa. O sea, la potencia de una prueba es la probabilidad de encontrar un resultado positivo dado que este existe. Una de las formas de representar la potencia de un test es a través del siguiente gráfico:
Del gráfico, es posible visualizar que a medida que aumenta la diferencia en la media de la población, se obtienen mayores valores de potencia estadística.
Recordada que es la potencia de un test de hipótesis, a continuación, usted deberá programar una función que sea capaz de obtener la potencia de un Z-test one-sample y two-sample. Para esto por favor considere los siguientes puntos:
- Crear una función que posea los siguientes argumentos:
function(n1=NULL, sigma1=0.5,
n2=NULL,sigma2=0.5, mu.Ha=0 ,
mu.True=0, alfa=0.05)De los argumentos, tendremos que: \(n1\) representa la cantidad de datos para la muestra 1, \(sigma1\) es la desviación estándar de la muestra 1, \(n2\) la cantidad de datos para la muestra 2, \(sigma2\) la desviación estándar para la muestra 2, \(mu.Ha\) el mu del test de hipótesis y \(mu.True\) la media de la población real. Notar que la presencia de una segunda muestra solo es para el caso two-sample, para el caso one-sample el argumento de entrada \(n2\) debería ser nulo.
- La función creada debe ser capaz de calcular el Z-test con el método One-sided (utilice solo la cola superior de la alternativa one-sided). Notar que la función al recibir un argumento nulo en \(n2\) debería asumir que se trata de un test one-sample automáticamente.
- Al recibir un valor no nulo para \(n2\), \(mu.Ha\) representará la diferencia entre las medias de las muestras y \(mu.True\) la diferencia de las medias de la población de las muestras 1 y 2.
- La salida de la función deberá retornar la potencia del test y un plot de las gaussianas que conforman el test de hipótesis. Para el caso del plot, observe los ejemplos de plot dispuestos más abajo.
- Si utiliza el esqueleto propuesto, complete y comente que realiza cada una de las partes de la función one-sample entregada.
Codificada la función realice los siguientes experimentos:
- Obtener el gráfico de potencia al variar la media poblacional para los siguientes argumentos de entrada:
\[ n1=16, sigma1=16, mu.Ha=100 , mu.True=Variar, alfa=0.05 \] \[ n1=16, sigma1=16, mu.Ha=100 , mu.True= Variar, alfa=0.01 \] \[ n1=16, sigma1=16, mu.Ha=100 , mu.True= Variar, alfa=0.1 \]
Se le recomienda que la variación se realice a través de un for y grafique las curvas dentro de un mismo gráfico para observar potenciales diferencias entre ellas.
- Diseñe un experimento Two-sample y visualice cómo se comportan las distribuciones normales de la hipótesis nula y la hipótesis alternativa al variar \(\alpha\).
Para el diseño de experimentos y/o comprobación de sus métodos puede serles útiles (no hay problema si decide utilizar los mismos ejemplos):
- one-sample: Power Functions
- Two-Sample: Simple Power Calculation for Two-Sample Z Test
Respuesta
# Power Function, El esqueleto posee como ejemplo como obtener la potencia de un z-test one-sample.
# Si utiliza este esqueleto deberá comentar la función que cumple cada una de las partes entregadas
power.z.test <- function(n1=NULL, sigma1=0.5,
n2=NULL, sigma2=0.5,
mu.Ha=0, mu.True=0, alfa=0.05){
if(is.null(n2)){
# Cálculo de power dado un sample
Z = qnorm(1-alfa)
denominador = sigma1/sqrt(n1)
X_bar = Z*denominador + mu.Ha
numerador = X_bar - mu.True
Z = numerador/denominador
Power = 1 - pnorm(Z)
# Se calculan los límites gráficos de power,
min_lim = min(rnorm(1000, mean=mu.Ha, sd=denominador)) -
round(min(rnorm(1000, mean=mu.Ha, sd=denominador)))%%10
max_lim = max(rnorm(1000, mean=mu.True, sd=denominador)) +
round(max(rnorm(1000, mean=mu.True, sd=denominador)))%%10
# Ploteo del power del test de hipotesis.
# En rojo queda la distribución del test de hipótesis, mientras que en azul la distribución real.
# Además, el área llenada en rojo es la visualización del power.
plot <- ggplot(data.frame(x = c(min_lim, max_lim)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = mu.Ha, sd = denominador),
col='red') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denominador),
col='blue') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denominador),
xlim = c(X_bar,max_lim), geom = "area", fill='red') +
geom_vline(xintercept = X_bar, linetype="dotted", size=1) +
annotate(x=X_bar, y=+Inf,label="alpha", vjust=2, geom="label") +
theme_minimal() +
ggtitle("H0 vs Ha") +
xlab(expression(bar(X))) + ylab("Density")
}
if(!is.null(n2)){
# Cálculo de power dado un sample
Z = qnorm(1-alfa)
denominador1 = sigma1/sqrt(n1)
denominador2 = sigma2/sqrt(n2)
denominador = sqrt(denominador2^2 + denominador1^2)
X_bar = Z*(denominador)
Z = (X_bar - mu.Ha)/denominador
Power = 1 - pnorm(Z)
# Se calculan los límites gráficos de power,
min_lim = min(rnorm(1000, mean=0, sd=denominador)) -
round(min(rnorm(1000, mean=0, sd=denominador)))%%10
max_lim = max(rnorm(1000, mean=mu.Ha, sd=denominador)) +
round(max(rnorm(1000, mean=mu.Ha, sd=denominador)))%%10
# Ploteo del power del test de hipotesis.
# En rojo queda la distribución del test de hipótesis, mientras que en azul la distribución 2
# Además, el área llenada en rojo es la visualización del power.
plot <- ggplot(data.frame(x = c(min_lim, max_lim)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = 0, sd = denominador),
col='red') +
stat_function(fun = dnorm, args = list(mean = mu.Ha, sd = denominador),
col='blue') +
stat_function(fun = dnorm, args = list(mean = mu.Ha, sd = denominador),
xlim = c(X_bar, max_lim), geom = "area", fill='red') +
geom_vline(xintercept = X_bar, linetype="dotted", size=1) +
annotate(x=X_bar, y=+Inf,label="alpha", vjust=2, geom="label") +
theme_minimal() +
ggtitle(paste("H0 vs Ha, with alpha = ", alfa, sep = "")) +
xlab("Means Difference") + ylab("Density")
}
# Como R no permite retornar dos salidas usamos una lista
# Los resultados se llaman con $plot o $power
return(list(plot=plot,power=Power))
}Experimento 1:
Para esto crearemos una función que genere los datos que necesitamos para los praficos.
# Función que entrega un dataframe de los power calculados para distintos
# valores del mu real. Este calcula el power para 20 valores de mu que van
# desde 0 a 2 desviaciones estándar del test de hipótesis
powers_realmu <- function(n1=NULL, sigma1=0.5,
mu.Ha=0, mu.True=0, alfa=0.05){
# Creo el dataframe con z y su power respectivo
plots <- data.frame(z = double(),
power = double())
# Creo un vector de los valores del experimento y de los valores de x
experiments = c()
x_vals = c()
step = sigma1*2/20
# Por cada step
for (i in 0:20) {
# Calculo el valor de la media real según el step
val = mu.Ha + step*i
experiments <- append(experiments, val)
x_vals = append(x_vals, i/10)
}
# Por cada valor de mu
for (value in experiments){
# Calculo el power utilizando los valores de mu
p <- power.z.test(n1, sigma1, NULL, 0, mu.Ha, value, alfa)
z <- data.frame(z = value, powers = p[2])
# Lo guardo en el dataframe
plots <- rbind(plots, z)
}
# Guardo tambien los valores de los steps
plots$x = x_vals
return(plots)
}Ahora basta realizar un for para los alphas sugeridos, guardar estos datos dentro del mismo dataframe y finalmente plotearlos de a cuerdo a su alpha como se muestra a continuación:
# Creo un dataframe para guardar los experimentos
vals <- data.frame(x = double(),
z = double(),
power = double(),
alpha = factor())
# Por cada valor de alpha, utilizo la función para variar el mu real
for (i in c(0.05, 0.01, 0.1)){
sample <- powers_realmu(n1=16, sigma1=16, mu.Ha=100, mu.True=value, alfa=i)
sample$alpha = factor(i)
vals <- rbind(vals, sample)
}
plot <- ggplot(vals, aes(x=z, y=power, color=alpha)) +
geom_line() +
geom_point() +
ylim(0, 1) +
xlab("Real mu values") +
ylab("Power") +
ggtitle("Power of the t Test \n n = 16, sigma = 16, mu = 100") +
theme(plot.title = element_text(hjust = 0.5))
plot
Notemos que al variar el valor de alpha las tres curvas comienzan en tres valores distintos. Esto se debe a que cuando la media real es idéntica a la del test de hipótesis estas se solapan perfectamente, por lo que el valor inicial cuando Real.mu = Ha.mu es el valor probabilístico de alpha en una distribución normal (1% en verde, 5% en rojo y 10% en azul).
Además de lo anterior es fácil ver que las tres curvas son identicas, pero trasladadas en el eje x dependiendo del valor inicial. Al “partir” de un alpha mayor, la probabilidad acumulada de rechazar el test de hipótesis aumenta más rápido que en el caso contrario, lo cuál se debe a la forma de la distribución normal y al grado de aumento de esta a medida que se mueve el mu real. En otras palabras, estamos en momentos distintos de la curva dependiendo del alpha.
Experimento 2:
Dado que ahora estamos en el caso two sample, podemos utilizar la misma función power.z.test para poder conseguir los datos necesarios:
p1 <- power.z.test(n1=50, sigma1=20, n2=50, sigma2=20, mu.Ha=10,
mu.True=10, alfa=0.025)[[1]]
p2 <- power.z.test(n1=50, sigma1=20, n2=50, sigma2=20, mu.Ha=10,
mu.True=10, alfa=0.05)[[1]]
p3 <- power.z.test(n1=50, sigma1=20, n2=50, sigma2=20, mu.Ha=10,
mu.True=10, alfa=0.1)[[1]]
p4 <- power.z.test(n1=50, sigma1=20, n2=50, sigma2=20, mu.Ha=10,
mu.True=10, alfa=0.2)[[1]]
grid.arrange(p1, p2, p3, p4,
# Numero de filas y columnas en los que los queremos sortear
ncol=2, nrow=2,
# Títulos
top = textGrob("Two Sample test with different alpha values \n n1=50, sigma1=20, n2=50, sigma2=20, \n mu1-mu2=10, x1-x2=10",
gp=gpar(fontsize=15)))
Claramente al variar el alfa cambiamos también la restricción de rechazar la hiótesis nula. A menor alpha, menor es la probabilidad de rechazar la hipótesis nula y por ende también es menor su power, y por otro lado al aumentar el valor de alpha este acerca el límite del alfa y aumenta la probabilidad de rechazar la hipótesis nula, lo cuál aumenta su power.
Pregunta 2: Z-test
Esta pregunta tiene como objetivo comprender como funciona un test de hipótesis y como deberíamos abordar la realización de múltiples test de hipótesis con datos reales.
La pregunta deberá ser desarrollada utilizando el dataset marketing_campaign.csv. Con esto, deberá programar un Z-test, con el cual estudiará a través de experimentos el Income de personas con los grados académicos Graduation, Master y PhD. Para realizar esto considere la elaboración de los siguientes puntos de forma secuencial:
- Modificar el dataframe entregado generando un estructura apta para el test de hipótesis. Una estructura que se les aconseja utilizar son vectores con los valores que representan a los grados académicos
Graduation,MasteryPhDpor separado.
Ejemplo de estructura
Por ejemplo para el caso de Graduation pueden generar estructuras de la siguiente forma:
| ID | Graduation |
|---|---|
| 5524 | 58138 |
| 2174 | 46344 |
| 4141 | 71613 |
| 6182 | 26646 |
| 965 | 55635 |
| … | … |
Donde los valores en la fila de Graduation representan los sueldos de las diferentes personas que conforman el dataset. Un punto importante a considerar es que los datos para los diferentes grados académicos poseen diferentes numero de datos (no se asusten por esto).
Programar el método Z-test con la metodología one sample y two sample, obteniendo los p-valores a través de las alternativas one-sided y two-sided. Para el caso de one-sided, cree una función capaz de obtener la cola menor y mayor de la gaussiana.
El calculo de las diferentes alternativas para calcular los p-valores deberá ser un argumento de su función, donde señalando ‘menor’,‘mayor’ (para los casos one-sided) y ‘two-sided’ deberá obtener el valor pertinente para cada caso.
Genere una función que permita realizar solo múltiples test del tipo two-sample y aplique bonferroni correction a los p-valores obtenidos. Notar que los múltiples test deberá realizar la comparación entre todos los elementos de entrada, por ejemplo si deseamos comparar los ingresos de
Graduation,MasteryPhD, se deberían comparar los ingresos deGraduationv/sMaster,Graduationv/sPhDyMasteryPhD
Codificada las funciones, realice los siguientes experimentos con su función de test de hipótesis:
Compruebe si la media de los ingresos para la variable
Graduationes similar a 52000. Señale formalmente este experimento y obtenga los p-valores para las alternativas one-sided y two-sided.Compruebe si la diferencia entre los ingresos de las personas con el grado académico
Graduationes cercana a cero en relación a la recibida por losMasteryPhD. Para este punto utilice la función que le permite realizar múltiples test del tipo two-sample.
Para los diferentes experimentos considere que la desviación estandar de la población para los diferentes income son los siguientes:
\[\sigma_{Graduation} = 28180\] \[\sigma_{Master} = 20160\] \[\sigma_{PhD} = 20615\]
Respuesta:
df = read.csv('marketing_campaign.csv', sep='\t')
# Implementación de Z-test one-sided y two-sided
# Puede utilizar este esqueleto
z_test <- function(data1=NULL, sigma1=0.5, data2=NULL, sigma2=0.5,
mu.Ha=0, test.type = c('one-sided','two-sided'),
verbose=TRUE){
if(length(test.type)>=2){
print("Por favor escoge un tipo de Test: ´one-sided´ o ´two-sided´ ")
return()
}
else if(length(test.type)==1 && !(test.type %in% c('menor','mayor','two-sided'))){
print("Por favor escoge un tipo de Test: ´menor´, ´mayor´ o ´two-sided´")
return()
}
else if(is.null(data2)){
# P-value
if(test.type=='menor'){
}
else if(test.type=='mayor'){
}
else if(test.type=='two-sided'){
}
# Texto de Salida
if(verbose){
cat("\tOne-sample Z-Test:\n\nData analizada:",
deparse(substitute(data1)), "\nZ=", Z_score,
"P-value=", p_value, "\n\n",sep=" ")
}
return(p_value)
}
else if(!is.null(data2)){
# Hypothesis test
# p-value
if(test.type=='menor'){
p_value = ...
}
else if(test.type=='mayor'){
p_value = ...
}
else if(test.type=='two-sided'){
p_value = ...
}
# Texto de Salida
if(verbose){
cat("\tTwo-sample Z-Test:\n\nData analizada:",
deparse(substitute(data1)),"y",
deparse(substitute(data2)), "\nZ=",
Z_score, "P-value=", p_value, "\n\n",sep=" ")
}
return(p_value)
}
}z.test.multiple_testing <- function(){
}Pregunta 3: Testeando multiples hipotesis y Bonferroni Correction
El objetivo de este problema es estudiar como realizar múltiples test de hipótesis simultáneamente. Para esto en primer lugar se estudiara el método “intuitivo”, donde veremos sus limitantes y se comparará con el método llamado Bonferroni correction, posteriormente se realizará un estudio practico con el dataset ratones.csv.
Un investigador se ha colocado en contacto con ustedes señalándoles que realiza diariamente test de hipótesis entre las muestras que toma día a día en su laboratorio. Con esto, al investigador le urge saber si realizar multiples test de hipótesis sin una corrección podría afectar la toma de decisiones. Para comprobar esto, les solicita comprobar matemáticamente como se comporta la probabilidad de obtener al menos un resultado significativos al azar de sus experimentos diarios. Para esto, les señala que la la probabilidad de obtener un experimento por azar puede ser simulado a través de los casos exitosos de una binomial (valores mayores a cero), donde el numero de observaciones son la cantidad de experimentos (\(m\)) y la probabilidad queda dada por \(\alpha\) del test.
A continuación, se entregan unas indicaciones mas especificas para desarrollar la pregunta:
- Complete el código presentado a continuación que le permite calcular la probabilidad empírica de que obtenga al menos un resultado significativo para significancia \(\alpha\) y cantidad de experimentos \(m\) arbitrarios.
- Se puede verificar que para un nivel de significancia \(\alpha\) y \(m\) experimentos independientes la probabilidad de que se tenga al menos un resultado significativo por azar es \[\mathbb{P}(\text{obtener al menos resultado significativo por azar})=1-(1-\alpha)^{m}\]
- Considere \(\alpha = 0.05\), grafique la probabilidad empírica y real variando el valor de \(m\) ¿Se parecen sus resultados? ¿Que sucede cuando la cantidad de experimentos crece mucho? ¿Este comportamiento depende del valor de significancia \(\alpha\)? ¿Es útil este método para la realización de múltiples test de hipótesis?
- Para solucionar los inconvenientes del método anterior es posible utilizar el método de Bonferroni correction, modifique su código anterior para verificar lo anterior ¿Mejoran los resultados? ¿cual podría ser un problema si es que \(m\) es muy grande?
- Ejecute el siguiente código que calcula el \(p\)-valor usual y el \(p\)-valor asociado a Bonferroni (que corresponde al \(p\)-valor * m donde \(m\) es el numero de experimentos), ¿Cuantos valores que originalmente se hubieran aceptado fueron rechazados si \(\alpha = 0.05\)? ¿Que implica esto sobre el nivel de falsos negativos de este metodo?
data <- read.csv("ratones.csv",sep= ";", stringsAsFactors = T)
head(data)Respuesta Aquí:
probEmpirica <- function(alpha, m){
n <- 1000 # Cantidad de veces que se va a repetir el experimento
res <- rbinom(n, m, alpha)
prob <- length(which(res>0)) / n
return(prob)
}Con nuestra función lista, podemos pasar a generar los datos necesarios para distintos valores de m:
probs <- data.frame(
m = integer(),
prob = double(),
type = integer())
prob <- 0.05
m <- c(1:100)
# Por cada tamaño de muestras
for(i in m) {
e <- probEmpirica(prob, i)
t <- 1-(1-prob)^i
z1 <- data.frame(m=i, prob=e, type=factor("Empirical"))
z2 <- data.frame(m=i, prob=t, type=factor("Theorical"))
probs <- rbind(probs, z1)
probs <- rbind(probs, z2)
}Correccion de Bonferroni:
Aquí solo debemos modificar la probabilidad de cada iteración de alpha a alpha/m con m el número de iteraciones.
prob <- 0.05
m <- c(1:100)
# Por cada tamaño de muestras
for(i in m) {
e <- probEmpirica(prob/i, i)
t <- 1-(1-prob/i)^i
z1 <- data.frame(m=i, prob=e, type=factor("Empirical B."))
z2 <- data.frame(m=i, prob=t, type=factor("Theorical B."))
probs <- rbind(probs, z1)
probs <- rbind(probs, z2)
}Ya obtenidos todos los valores, utilizando las probabilidades con y sin la corrección de bonferroni, podemos graficarlos como sigue:
# Generamos el gráfico con los resultados
g <- ggplot(probs, aes(m, prob)) +
geom_line(aes(colour=type)) +
ylim(0, 1) +
xlim(0, 100) +
xlab("m") +
ylab("Probabilidad") +
ggtitle("Probabilidad teorica y empirica de multiples tests\ncon y sin la corrección de Bonferroni\nalpha = 0.05") +
theme(plot.title = element_text(hjust = 0.5))
g
Claramente los valores epmíricos y teóricos son muy similares (donde los empíricos presentan ruido dentro de un margen de error). Sin mebargo, entre los valores con y sin la corrección de Bonferroni tenemos una clara diferencia de valores, donde los que tienen la corrección se mantienen constantes al rededor del valor 0.05, y los que no tienen la corrección son crecientes y alcanzan su máximo con un valor aprocimado de m=80.
Lo anterior nos indica que a mayor cantidad de tests múltiples, mayor es el error probabilístico contenido, y que este crece muy rápidamente a medida que la cantidad de tests aumenta. Con la correción de Bonferroni podemos asegurar que los errores de este tipo (falsos positivos) se mantengan estables al rededor del valor de alpha.
Claramente la forma de la probabilidad acumulada depende del valor de alpha, ya que a mayor alpha más rápido crecerá la curva sin la corrección de Bonferroni, y a su vez, el valor al que se estabiliza la probabilidad utilizando la corrección es mayor.
En cuanto a si este método es útil, la respuesta es que depende del tipo de test de hipótesis. Si queremos que el error no aumente cuando estamos trabajando con errores de tipo I, entonces este método sí funciona ya que mantenemos constante la probabilidad de valores significativos, y por ende, tenemos menor error. Sin embargo, si queremos disminuir la probabilidad de valores significativos en caso de errores de tipo II (falsos negativos), entonces esta corrección es inválida, ya que tiene el efecto inverso.
Pregunta 4: Regression Lineal sin comandos.
El objetivo de la siguiente pregunta es aplicar los conceptos de regresión lineal vistos en clases para implementar desde cero un función capaz de realizar una regresión simple y múltiple.
Para este problema, ustedes deberán estudiar el comportamiento de los clientes de un holding de salud. Para esto, se les hace entrega del dataset insurance.csv para que estudien la creación de un modelo lineal con sus datos. Antes de comenzar a trabajar, se señalan las diferentes variables que componen al dataset:
- age: Señala la edad de cada uno de los sujetos.
- sex: Si es mujer es igual a 1, si es hombre es igual a 0.
- bmi: Indice de masa corporal del cliente.
- children: Señala cuantos hijos tiene cada uno de los sujetos.
- smoker: Variable binaria que cuando es 1 señala que el cliente es fumador (0 en caso contrario).
- charges: Gastos médicos de cada uno de los clientes.
Es importante que considere que cada una de las filas representa un cliente distinto para el holding.
Dentro del estudio, el holding de salud le solicita estudiar los comportamientos de los clientes fumadores y no fumadores, por lo que se le aconseja separar el dataframe original en fumadores y no fumadores. En el estudio, realicen un modelo lineal que tiene como variable de respuesta a charges y los datos que mejor se correlacionan para los clientes fumadores y no fumadores. Para esto, deberán realizar las siguientes actividades:
Parte I
- Programe un modelo lineal simple escogiendo la variable numérica que tiene mayor relación con la variable de respuesta. Recuerde justificar la elección de la variable numérica cuantitativamente.
- Señale tanto el \(R^2\) como el \(R^2-adjustado\) del modelo.
- Grafique el scatterplot de los datos y la linea que ajusta a la regresión lineal obtenida.
Parte II
- Entrene un modelo lineal multivariable escogiendo dos variables numéricas que posean la mayor relación con
charges. - Estudie si el modelo multivariable posee mejor desempeño que el modelo simple y comente los resultados. ¿Es recomendable la utilización de los modelos creados para la predicción de nuevas entradas?. Para este análisis puede utilizar los valores de test de hipótesis entregados por el comando
lm(), ya que esto le servirá para observar si la regresión lineal es significativa.
Nota: No esta permitido utilizar comandos que obtengan los valores solicitados directamente a menos que se le permita en la pregunta.
Respuesta
Parte I
Primero, leemos la tabla
data <- read.csv("insurance.csv", stringsAsFactors = T)
head(data)Separaremos el dataframe en dos, donde uno representa los fumadores y el otro los no fumadore, y calculamos la correlación entre las variables por cada grupo para determinar la variable que mejor se relacione con charges.
# Calculo de la matriz de correlacion para su visualización
smokers <- data[data$smoker == "yes",]
non_smokers <- data[data$smoker == "no",]
correlation_smokers = cor(smokers[c(1,3,4,7)])
correlation_non_smokers = cor(non_smokers[c(1,3,4,7)])
print("smokers:")## [1] "smokers:"round(correlation_smokers, 3)## age bmi children charges
## age 1.000 0.060 0.081 0.368
## bmi 0.060 1.000 -0.013 0.806
## children 0.081 -0.013 1.000 0.036
## charges 0.368 0.806 0.036 1.000print("non smokers:")## [1] "non smokers:"round(correlation_non_smokers, 3)## age bmi children charges
## age 1.000 0.123 0.033 0.628
## bmi 0.123 1.000 0.019 0.084
## children 0.033 0.019 1.000 0.139
## charges 0.628 0.084 0.139 1.000Podemos ver claramente que dependiendo de si se es fumador o no, la variable numérica que mejor representa a cada grupo es distinta. En el caso de los fumadores la mejor métrica es el bmi, mientras que en el caso de los no fumadores es la edad.
Ahora, creamos una función que calcule la regresión lineal simple dados dos vectores. Con X la variable independiente e Y la variable de respuesta:
lm_reg <- function(X, Y){
X_bar = mean(X)
Y_bar = mean(Y)
num = 0
den = 0
for (i in 1:length(X)){
num = num + (X[i]-X_bar)*(Y[i]-Y_bar)
den = den + (X[i]-X_bar)^2
}
b1 = num/den
b0 = Y_bar - b1*X_bar
y_gorro = c()
for (i in 1:length(X)){
y_gorro = c(y_gorro, (X[i]*b1 + b0))
}
R2 = cor(Y, y_gorro)
adjusted_R2 = (1 - ((1-R2)*(length(X)-1))/(length(X)-2))
return(list(b0, b1, R2, adjusted_R2))
}Ahora, crearemos dos listas de resultados con ambos datasets, el de fumadores y el de no fumadores
lm_smokers <- lm_reg(smokers$bmi, smokers$charges)
lm_non_smokers <- lm_reg(non_smokers$age, non_smokers$charges)print(paste("smokers: R2=", lm_smokers[[3]], " adjusted_R2=", lm_smokers[[4]]))## [1] "smokers: R2= 0.806480607015541 adjusted_R2= 0.80576913865898"print(paste("non smokers: R2=", lm_non_smokers[[3]], " adjusted_R2=", lm_non_smokers[[4]]))## [1] "non smokers: R2= 0.62794678376642 adjusted_R2= 0.627596451171096"plot <- ggplot(smokers, aes(x=bmi, y=charges)) +
geom_point() +
geom_abline(intercept = lm_smokers[[1]], slope = lm_smokers[[2]], color="red") +
xlab("BMI") +
ylab("Charges") +
ggtitle("Simple Linear Regression for Smokers\nCharges vs BMI") +
theme(plot.title = element_text(hjust = 0.5))
plot
plot <- ggplot(non_smokers, aes(x=age, y=charges)) +
geom_point() +
geom_abline(intercept = lm_non_smokers[[1]], slope = lm_non_smokers[[2]], color="blue") +
xlab("Ages") +
ylab("Charges") +
ggtitle("Simple Linear Regression for Non Smokers\nCharges vs Age") +
theme(plot.title = element_text(hjust = 0.5))
plot
Parte II
Para esta parte basta crear una nueva funcion de regresión lineal pero con dos variables. Ahora X e Y son las variables independientes y la variable Z la variable de respuesta:
lm_reg2 <- function(X1, X2, Z){
X1_bar = mean(X1)
X2_bar = mean(X2)
Z_bar = mean(Z)
sum_x22 = 0
sum_x12 = 0
sum_x1x2 = 0
sum_x1z = 0
sum_x2z = 0
for (i in 1:length(X1)){
sum_x22 = sum_x22 + (X2[i]-X2_bar)^2
sum_x12 = sum_x12 + (X1[i]-X1_bar)^2
sum_x1x2 = sum_x1x2 + (X1[i]-X1_bar)*(X2[i]-X2_bar)
sum_x1z = sum_x1z + (X1[i]-X1_bar)*(Z[i]-Z_bar)
sum_x2z = sum_x2z + (X2[i]-X2_bar)*(Z[i]-Z_bar)
}
b1 = (sum_x22*sum_x1z - sum_x1x2*sum_x2z)/(sum_x12*sum_x22 - sum_x1x2^2)
b2 = (sum_x12*sum_x2z - sum_x1x2*sum_x1z)/(sum_x12*sum_x22 - sum_x1x2^2)
a = Z_bar - b1*X1_bar - b2*X2_bar
z_gorro = c()
for (i in 1:length(X1)){
z_gorro = c(z_gorro, (a + b1*X1[i] + b2*X2[i]))
}
R2 = cor(Z, z_gorro)
adjusted_R2 = (1 - ((1-R2)*(length(X1)-1))/(length(X1)-3))
return(list(a, b1, b2, R2, adjusted_R2))
}Ahora calculando nuevamente la regresión para los casos de fumador y no fumador, pero con las dos correlaciones más grandes de cada uno tendremos:
lm_smokers2 <- lm_reg2(smokers$bmi, smokers$age, smokers$charges)
lm_non_smokers2 <- lm_reg2(non_smokers$age, non_smokers$children, non_smokers$charges)print(paste("smokers: R2=", lm_smokers2[[4]], " adjusted_R2=", lm_smokers2[[5]]))## [1] "smokers: R2= 0.867894155513266 adjusted_R2= 0.866919204631445"print(paste("non smokers: R2=", lm_non_smokers2[[4]], " adjusted_R2=", lm_non_smokers2[[5]]))## [1] "non smokers: R2= 0.638941967880991 adjusted_R2= 0.63826136838595"Ya viendo únicamente estos valores de R2 y R2 ajustados tenemos que la regresión efectivamente tiene mejor desempeño que utilizando una única variable independiente. Esta mejora es más significativa en el caso de la tabla de fumadores (6-7 puntos de mejora), mientras que para los no-fumadores esta mejora en el desempeño no es tan grande (1-2 puntos). Esto podría indicar que la edad y la variable añadida “children” no aporta mayor información en el caso de los datos de personas no fumadoras, mientras que en el caso de los fumadores podemos concluir que la variable bmi y la edad si pueden describir de buena manera el comportamiento de los datos.
linear_s <- lm(smokers$charges ~ smokers$bmi + smokers$age)
linear_ns <- lm(non_smokers$charges ~ non_smokers$age + non_smokers$children)
summary(linear_s)##
## Call:
## lm(formula = smokers$charges ~ smokers$bmi + smokers$age)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14604.4 -4315.1 -240.5 3638.0 29316.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22367.45 1931.86 -11.58 <2e-16 ***
## smokers$bmi 1438.09 55.22 26.05 <2e-16 ***
## smokers$age 266.29 25.06 10.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5754 on 271 degrees of freedom
## Multiple R-squared: 0.7532, Adjusted R-squared: 0.7514
## F-statistic: 413.6 on 2 and 271 DF, p-value: < 2.2e-16summary(linear_ns)##
## Call:
## lm(formula = non_smokers$charges ~ non_smokers$age + non_smokers$children)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2548.4 -1894.6 -1364.5 -680.8 24490.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2658.80 435.44 -6.106 1.43e-09 ***
## non_smokers$age 265.57 10.06 26.408 < 2e-16 ***
## non_smokers$children 581.06 116.27 4.998 6.79e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4615 on 1061 degrees of freedom
## Multiple R-squared: 0.4082, Adjusted R-squared: 0.4071
## F-statistic: 366 on 2 and 1061 DF, p-value: < 2.2e-16Ya comparandolo con el modelo generado por R, podemos ver que nuestra modelo si tiene un comportamiento similar y significativo. Además, podemos utilizar las herramientas que nos entrega R para modelar los datos de cada dataset, y la regresión calculada en forma de plano:
### Calculate z on a grid of x-y values
x1.seq <- seq(min(smokers$bmi),max(smokers$bmi),length.out=25)
x2.seq <- seq(min(smokers$age),max(smokers$age),length.out=25)
z <- t(outer(x1.seq, x2.seq, function(x,y) lm_smokers2[[1]] + lm_smokers2[[2]]*x + lm_smokers2[[3]]*y))
#### Draw the plane with "plot_ly" and add points with "add_trace"
p <- plot_ly(x=~x1.seq, y=~x2.seq, z=~z, colors = c("#f5cb11", "#cf2525"), opacity=0.5, type="surface") %>%
add_trace(x=smokers$bmi, y=smokers$age, z=smokers$charges, mode="markers", type="scatter3d",
marker = list(opacity=0.7, symbol=105, color="#f03c3c")) %>%
layout(title = 'Regression Plane for smokers', scene = list(
aspectmode = "manual",
xaxis = list(title = "bmi"),
yaxis = list(title = "age"),
zaxis = list(title = "charge")))
p### Calculate z on a grid of x-y values
x1.seq <- seq(min(non_smokers$age),max(non_smokers$age),length.out=25)
x2.seq <- seq(min(non_smokers$children),max(non_smokers$children),length.out=25)
z <- t(outer(x1.seq, x2.seq, function(x,y) lm_non_smokers[[1]] + lm_non_smokers[[2]]*x + lm_non_smokers[[3]]*y))
#### Draw the plane with "plot_ly" and add points with "add_trace"
p2 <- plot_ly(x=~x1.seq, y=~x2.seq, z=~z, colors = c("#a6ffb0", "#4237db"), opacity=0.5, type="surface") %>%
add_trace(x=non_smokers$age, y=non_smokers$children, z=non_smokers$charges, mode="markers", type="scatter3d",
marker = list(opacity=0.7, symbol=105, color="#110ec2")) %>%
layout(title = 'Regression Plane for non smokers', scene = list(
aspectmode = "manual",
xaxis = list(title = "age"),
yaxis = list(title = "children"),
zaxis = list(title = "charge")))
p2A work by CC6104